InitDischargeRouting Subroutine

public subroutine InitDischargeRouting(fileini, time, path_out, flowdirection, domain, storage, netRainfall, dtrk)

Initialize surface routing

Arguments

Type IntentOptional Attributes Name
character(len=*), intent(in) :: fileini

name of configuration file
type(DateTime), intent(in) :: time
character(len=*), intent(in) :: path_out

path of output folder
type(grid_integer), intent(in) :: flowdirection
type(grid_integer), intent(in) :: domain
type(grid_real), intent(out) :: storage
type(grid_real), intent(out) :: netRainfall
integer(kind=short), intent(out) :: dtrk

Variables

Type Visibility Attributes Name Initial
real(kind=float), public :: area

basin area (m2)
character(len=30), public :: channelInitiationMethod

method to assign channel initiation
real(kind=float), public :: channelInitiationThreshold

channel initiation threshold (m2)
real(kind=float), public :: channelValue
type(Diversion), public, POINTER :: d

point to current diversion channel
integer(kind=short), public :: exportChannelGrid

set channel grid export
real(kind=float), public :: hillslopeBankSlope

hillslope bank slope (degree)
real(kind=float), public :: hillslopeKs

hillslope roughness (m^1/3 s-1)
real(kind=float), public :: hillslopelWidth

hillslope section width (m)
integer(kind=short), public :: i
type(IniList), public :: ini
integer(kind=short), public :: j
integer(kind=short), public :: k
type(grid_integer), public :: maskTmp
integer(kind=short), public :: nmasks

number of masks to assign routing parameters
type(Reservoir), public, POINTER :: p

point to current reservoir
real(kind=float), public :: scalar
character(len=100), public :: stringTmp

Source Code

SUBROUTINE InitDischargeRouting  &
!
(fileini, time, path_out, flowdirection, domain, & 
  storage, netRainfall, dtrk )

IMPLICIT NONE

!Arguments with intent in:
CHARACTER (LEN = *), INTENT(IN) :: fileini !!name of configuration file
TYPE(DateTime), INTENT(IN)      :: time
CHARACTER (LEN = *), INTENT(IN) :: path_out  !!path of output folder
TYPE(grid_integer), INTENT(IN)  :: flowdirection
TYPE(grid_integer), INTENT(IN)  :: domain

!Arguments with intent out:
TYPE(grid_real), INTENT(OUT)          :: storage  !water stored in channel cell
TYPE(grid_real), INTENT(OUT)          :: netRainfall ! net rainfall rate [m/s]
INTEGER(KIND = short), INTENT(OUT)    :: dtrk

!local declarations:
TYPE (IniList) :: ini
INTEGER (KIND = short) :: i, j, k 
!CHARACTER (LEN=5)     :: ic !!option for initial condition
TYPE(Reservoir), POINTER :: p !!point to current reservoir
TYPE(Diversion), POINTER :: d !!point to current diversion channel
CHARACTER (LEN = 30) :: channelInitiationMethod !!method to assign channel initiation
REAL (KIND = float) :: channelInitiationThreshold !!channel initiation threshold (m2)
REAL (KIND = float) :: hillslopelWidth !!hillslope section width (m)
REAL (KIND = float) :: hillslopeKs !!hillslope roughness (m^1/3 s-1)
REAL (KIND = float) :: hillslopeBankSlope !!hillslope bank slope (degree)
REAL (KIND = float) :: area !! basin area (m2)
REAL (KIND = float) :: channelValue
INTEGER (KIND = short) :: exportChannelGrid !!set channel grid export
REAL (KIND = float) :: scalar
INTEGER (KIND = short) :: nmasks !! number of masks to assign routing parameters
TYPE (grid_integer) :: maskTmp
CHARACTER (LEN = 100) :: stringTmp

!-------------------------------end of declaration-----------------------------

!check stream network

IF ( .NOT. streamNetworkCreated ) THEN
   CALL Catch ('error', 'DischargeRouting',   &
          'stream network is not available, check morphological properties' )
END IF

!load configuration file
CALL IniOpen (fileini, ini)

!number of masks to assign discharge routing parameters
IF ( KeyIsPresent ('masks-number', ini) ) THEN
    nmasks = IniReadInt ('masks-number', ini)
ELSE
    CALL Catch ('error', 'DischargeRouting',   &
			   'masks-number missing in configuration file' )
END IF


!channel initiation
!------------------

!allocate channel grid
CALL NewGrid ( channel, domain, 0)

!set option to export channel grid
exportChannelGrid = IniReadInt ('export-channel-grid', ini)

!first use base mask

channelInitiationMethod = IniReadString ('channel-initiation-method', &
                                         ini, section = 'base-mask')
channelInitiationThreshold = IniReadReal ('channel-initiation-threshold', &
                                         ini, section = 'base-mask')

CALL ChannelInitiation (method = channelInitiationMethod, &
                           threshold = channelInitiationThreshold, &
                           mask = domain, flowAcc = flowAccumulation,&
                           flowDir = flowDirection, dem = dem, &
                           channel = channel ) 


! modify channel initiation grid if additional masks are required
DO k = 1, nmasks - 1
    
    stringTmp = 'mask' // ToString (k)
    
    channelInitiationMethod = IniReadString ('channel-initiation-method', &
                                         ini, section = stringTmp)
    channelInitiationThreshold = IniReadReal ('channel-initiation-threshold', &
                                         ini, section = stringTmp)
    
    CALL GridByIni (ini, maskTmp, section = stringTmp)

    CALL ChannelInitiation (method = channelInitiationMethod, &
                           threshold = channelInitiationThreshold, &
                           mask = maskTmp, flowAcc = flowAccumulation,&
                           flowDir = flowDirection, dem = dem, &
                           channel = channel ) 
    CALL GridDestroy (maskTmp)
    
END DO


!export channel grid
IF ( exportChannelGrid == 1 ) THEN
    CALL ExportGrid (channel, 'channel.asc', ESRI_ASCII)
END IF

!create discharge parameter maps
!-------------------------------

!initialize empty maps
CALL NewGrid (layer = manning, grid = domain, initial = 0. )
CALL NewGrid (layer = bankSlope, grid = domain, initial = 0. )
CALL NewGrid (layer = sectionWidth, grid = domain, initial = 0. )

!read channel parameter table(s)
CALL TableNew ( fileini, channelParameters )

!assign parameters to base mask
hillslopelWidth = IniReadReal ('hillslope-width', ini, section = 'base-mask')
hillslopeKs = IniReadReal ('hillslope-ks', ini, section = 'base-mask')
hillslopeBankSlope = IniReadReal ('hillslope-alpha', ini, section = 'base-mask')
 

DO i = 1, domain % idim
    DO j = 1, domain % jdim
        IF ( domain % mat (i,j) /= domain % nodata ) THEN
            area = flowAccumulation % mat (i,j) * CellArea (domain,i,j)
            IF ( channel % mat (i,j) == 0 ) THEN !hillslope cell
                manning % mat (i,j) = 1. / hillslopeKs
                bankSlope % mat (i,j) = hillslopeBankSlope
                sectionWidth % mat (i,j) = hillslopelWidth
            ELSE !channel cell
                
                !set manning roughness coefficient
                CALL TableGetValue ( valueIn =  area, &
                                     tables = channelParameters, &
                                     id = 'base-mask',  keyIn = 'threshold', &
                                     keyOut ='ks', match = 'linear', &
                                     valueOut = channelValue, &
                                     bound = 'extendlinear' )
                
                manning % mat (i,j) = 1. / channelValue
                
                !set river bank slope
                CALL TableGetValue ( valueIn =  area, &
                                     tables = channelParameters, &
                                     id = 'base-mask',  keyIn = 'threshold', &
                                     keyOut ='alpha', match = 'linear', &
                                     valueOut = channelValue, &
                                     bound = 'extendlinear' )
                
                bankSlope % mat (i,j) = channelValue
                
                !set section width
                CALL TableGetValue ( valueIn =  area, &
                                     tables = channelParameters, &
                                     id = 'base-mask',  keyIn = 'threshold', &
                                     keyOut ='width', match = 'linear', &
                                     valueOut = channelValue, &
                                     bound = 'extendlinear' )
                
                sectionWidth % mat (i,j) = channelValue
                
            END IF
        END IF
    END DO
END DO


!assign parameters to sub masks

DO k = 1, nmasks - 1
    
    stringTmp = 'mask' // ToString (k)

    hillslopelWidth = IniReadReal ('hillslope-width', ini, &
                                   section = stringTmp)
    hillslopeKs = IniReadReal ('hillslope-ks', ini, &
                               section = stringTmp)
    hillslopeBankSlope = IniReadReal ('hillslope-alpha', ini, &
                                       section = stringTmp)
    
    CALL GridByIni (ini, maskTmp, section = stringTmp)
 
    DO i = 1, domain % idim
        DO j = 1, domain % jdim
            IF ( maskTmp % mat (i,j) /= maskTmp % nodata ) THEN
                area = flowAccumulation % mat (i,j) * CellArea (domain,i,j)
                IF ( channel % mat (i,j) == 0 ) THEN !hillslope cell
                    manning % mat (i,j) = 1. / hillslopeKs
                    bankSlope % mat (i,j) = hillslopeBankSlope
                    sectionWidth % mat (i,j) = hillslopelWidth
                ELSE !channel cell
                
                    !set manning roughness coefficient
                    CALL TableGetValue ( valueIn =  area, &
                                         tables = channelParameters, &
                                         id = stringTmp,  keyIn = 'threshold', &
                                         keyOut ='ks', match = 'linear', &
                                         valueOut = channelValue, &
                                         bound = 'extendlinear' )
                
                    manning % mat (i,j) = 1. / channelValue
                
                    !set river bank slope
                    CALL TableGetValue ( valueIn =  area, &
                                         tables = channelParameters, &
                                         id = stringTmp,  keyIn = 'threshold', &
                                         keyOut ='alpha', match = 'linear', &
                                         valueOut = channelValue, &
                                         bound = 'extendlinear' )
                
                    bankSlope % mat (i,j) = channelValue
                
                    !set section width
                    CALL TableGetValue ( valueIn =  area, &
                                         tables = channelParameters, &
                                         id = stringTmp,  keyIn = 'threshold', &
                                         keyOut ='width', match = 'linear', &
                                         valueOut = channelValue, &
                                         bound = 'extendlinear' )
                
                    sectionWidth % mat (i,j) = channelValue
                
                END IF
            END IF
        END DO
    END DO
    
    
    CALL GridDestroy (maskTmp)
    
END DO

!initial condition
!-----------------

!input discharge
IF (SectionIsPresent('discharge-in', ini )) THEN
     IF (KeyIsPresent ('scalar', ini, 'discharge-in') ) THEN
        scalar = IniReadReal ('scalar', ini, 'discharge-in')
        CALL NewGrid (Pin, domain, scalar)
    ELSE
        CALL GridByIni (ini, Pin, section = 'discharge-in')
    END IF
ELSE !grid is optional: set to default = 0
   CALL NewGrid ( Pin, domain, 0. )
END IF

!output discharge
IF (SectionIsPresent('discharge-out', ini )) THEN
     IF (KeyIsPresent ('scalar', ini, 'discharge-out') ) THEN
        scalar = IniReadReal ('scalar', ini, 'discharge-out')
        CALL NewGrid (Pout, domain, scalar)
    ELSE
        CALL GridByIni (ini, Pout, section = 'discharge-out')
    END IF
ELSE !grid is optional: set to default = 0
   CALL NewGrid ( Pout, domain, 0. )
END IF

!lateral discharge
IF (SectionIsPresent('discharge-lat', ini )) THEN
     IF (KeyIsPresent ('scalar', ini, 'discharge-lat') ) THEN
        scalar = IniReadReal ('scalar', ini, 'discharge-lat')
        CALL NewGrid (Plat, domain, scalar)
    ELSE
        CALL GridByIni (ini, Plat, section = 'discharge-lat')
    END IF
ELSE !grid is optional: set to default = 0
   CALL NewGrid ( Plat, domain, 0. )
END IF

!allocate variables
CALL NewGrid (layer = Qin, grid = domain, initial = 0. )
CALL NewGrid (layer = Qout, grid = domain, initial = 0. )
CALL NewGrid (layer = storage, grid = domain, initial = 0. )
 
!net rainfall grid
CALL NewGrid (layer = netRainfall, grid = domain, initial = 0.)

!water depth grid
CALL NewGrid (layer = waterDepth, grid = domain, initial = 0.)

!topwidth 
CALL NewGrid (layer = topWidth, grid = domain, initial = 0.)

!diversions
IF (SectionIsPresent('diversions', ini)) THEN
    
    dtDiversion = IniReadInt ('dt', ini, section = 'diversions')
    IF ( dtDiversion /= 0 ) THEN
        !set dtDiversion = dtDischargeRouting
        dtDiversion = dtDischargeRouting
    END IF
    
    IF (dtDiversion > 0) THEN
        CALL InitDiversions ( IniReadString('file', ini, &
                              section = 'diversions') )
        
        !create grid with diversion location
        CALL NewGrid (divChannels, domain, domain % nodata ) 
        
        d => diversionChannels
        DO 
          i = d % r
          j = d % c
          divChannels % mat (i,j) = d % id
          
          IF ( channel % mat (i,j) == 0 ) THEN
             CALL Catch ('warning', 'DischargeRouting',   &
			   'diversion is not on channel: ', argument = ToString(d%id) )
          END IF
          
            
          IF ( .NOT. ASSOCIATED (d % next) ) THEN !found last element of list
              EXIT
          END IF
          
          !set next diversion
          d => d % next
        END DO
        
        !initialize files for output
        IF ( KeyIsPresent ('dt-out', ini, 'diversions' ) ) THEN
            
            dtOutDiversion = IniReadInt ('dt-out', ini, 'diversions')
            
            IF ( dtOutDiversion > 0) THEN
               timeDiversionsExport = time
               CALL OutDiversionsInit ( path_out )
            END IF
            
        END IF
    
    END IF
    
ELSE !section is mandatory: stop the program if missing
   CALL Catch ('error', 'DischargeRouting',   &
			   'error loading diversions: ' ,  &
			    argument = 'section not defined in ini file' )
END IF

!reservoirs
IF (SectionIsPresent('reservoirs', ini)) THEN
    
    !create grid with reservoirs location
    CALL NewGrid (dams, domain, domain % nodata ) 
    
    
    dtrk = IniReadInt ('dt', ini, section = 'reservoirs')
    IF (dtrk > 0) THEN
        CALL InitReservoirs (IniReadString('file', ini, &
                             section = 'reservoirs'), &
                             time, domain, pools)
        
        p => pools
        DO 
          i = p % r
          j = p % c
          dams % mat (i,j) = p % id
          
          IF ( channel % mat (i,j) == 0 ) THEN
             CALL Catch ('warning', 'DischargeRouting',   &
			   'reservoir is not on channel: ', argument = ToString(p%id) )
          END IF
          
            
          IF ( .NOT. ASSOCIATED (p % next) ) THEN !found last element of list
              EXIT
          END IF
  
           p => p % next
        END DO
    
            
        !initialize files for output
        IF ( KeyIsPresent ('dt-out', ini, 'reservoirs' ) ) THEN
            
            dtOutReservoirs = IniReadInt ('dt-out', ini, 'reservoirs')
            
            IF ( dtOutReservoirs > 0) THEN
               timePoolsExport = time
               CALL OutReservoirsInit (pools, path_out)
            END IF
            
        END IF
    
    END IF
    

ELSE !section is optional
    
   dtrk = 0
   
 END IF
 
!finished configuration, deallocate memory
CALL IniClose (ini)


END SUBROUTINE InitDischargeRouting